encode country,gen(id)
xtset id year
*一键取对数
foreach var of varlist ET RENT PR FR ER  FC TL IL {
gen ln`var' = ln(`var')
}

*描述性分析
tabstat ET RENT PR FR ER FC TL IL,stats(n mean median sd min max) c(s) f(%10.3f)
*导出描述性分析
logout,save(descripitive statistics) word replace:tabstat ET PR FR ER FC TL IL,stats(n mean median sd min max) c(s) f(%10.3f)
*LLC检验
xtunitroot llc lnET,trend demean
xtunitroot llc lnRENT,trend demean
xtunitroot llc lnPR,trend demean
xtunitroot llc lnFR,trend demean
xtunitroot llc lnER,trend demean
xtunitroot llc lnFC,trend demean
xtunitroot llc lnTL,trend demean
xtunitroot llc lnIL,trend demean
*一阶差分
gen dET=d.lnET
gen dRENT=d.lnRENT
gen dPR=d.lnPR
gen dFR=d.lnFR
gen dER=d.lnER
gen dFC=d.lnFC
gen dTL=d.lnTL
gen dIL=d.IL
xtunitroot llc d.lnET,trend demean
xtunitroot llc d.lnRENT,trend demean
xtunitroot llc d.lnPR,trend demean
xtunitroot llc d.lnFR,trend demean
xtunitroot llc d.lnER,trend demean
xtunitroot llc d.lnFC,trend demean
xtunitroot llc d.lnTL,trend demean
xtunitroot llc d.lnIL,trend demean
*IPS检验
xtunitroot ips lnET,trend demean
xtunitroot ips lnRENT,trend demean
xtunitroot ips lnPR,trend demean
xtunitroot ips lnFR,trend demean
xtunitroot ips lnER,trend demean
xtunitroot ips lnFC,trend demean
xtunitroot ips lnTL,trend demean
xtunitroot ips lnIL,trend demean

xtunitroot ips d.lnET,trend demean
xtunitroot ips d.lnRENT,trend demean
xtunitroot ips d.lnPR,trend demean
xtunitroot ips d.lnFR,trend demean
xtunitroot ips d.lnER,trend demean
xtunitroot ips d.lnFC,trend demean
xtunitroot ips d.lnTL,trend demean
xtunitroot ips d.lnIL,trend demean

*HT检验
xtunitroot ht lnET,trend demean
xtunitroot ht lnRENT,trend demean
xtunitroot ht lnPR,trend demean
xtunitroot ht lnFR,trend demean
xtunitroot ht lnER,trend demean
xtunitroot ht lnFC,trend demean
xtunitroot ht lnTL,trend demean
xtunitroot ht lnIL,trend demean
*一阶差分
xtunitroot ht d.lnET,trend demean
xtunitroot ht d.lnRENT,trend demean
xtunitroot ht d.lnPR,trend demean
xtunitroot ht d.lnFR,trend demean
xtunitroot ht d.lnER,trend demean
xtunitroot ht d.lnFC,trend demean
xtunitroot ht d.lnTL,trend demean
xtunitroot ht d.lnIL,trend demean
*协整检验
xtcointtest westerlund lnET lnRENT lnFR lnER lnPR  lnFC lnTL lnIL,trend

*门槛效应模型
*PR为解释变量
set seed 123
xthreg lnET lnER lnFR  lnFC lnTL lnIL,rx(lnPR) qx(lnRENT) thnum(1) bs(300) trim(0.01) grid(300) 
set seed 123
xthreg lnET lnER lnFR  lnFC lnTL lnIL,rx(lnPR) qx(lnRENT) thnum(2) bs(300 300) trim(0.01 0.01) grid(300) 
set seed 123
xthreg lnET lnER lnFR  lnFC lnTL lnIL,rx(lnPR) qx(lnRENT) thnum(3) bs(300 300 300) trim(0.01 0.01 0.05) grid(300) 
*FR为解释变量
set seed 123
xthreg lnET lnER lnPR  lnFC lnTL lnIL,rx(lnFR) qx(lnRENT) thnum(1) bs(300) trim(0.01) grid(300) 
set seed 123
xthreg lnET lnER lnPR lnFC lnTL lnIL,rx(lnFR) qx(lnRENT) thnum(2) bs(300 300) trim(0.01 0.01) grid(300) 
set seed 123
xthreg lnET lnER lnPR lnFC lnTL lnIL,rx(lnFR) qx(lnRENT) thnum(3) bs(300 300 300) trim(0.01 0.01 0.05) grid(300) 
*ER为解释变量
set seed 123
xthreg lnET lnPR lnFR lnFC lnTL lnIL,rx(lnER) qx(lnRENT) thnum(1) bs(300) trim(0.01) grid(300) 
set seed 123
xthreg lnET lnPR lnFR lnFC lnTL lnIL,rx(lnER) qx(lnRENT) thnum(2) bs(300 300) trim(0.01 0.01) grid(300) 
set seed 123
xthreg lnET lnPR lnFR lnFC lnTL lnIL,rx(lnER) qx(lnRENT) thnum(3) bs(300 300 300) trim(0.01 0.01 0.05) grid(300) 

**绘制双门槛图:
_matplot e(LR21), columns(1 2) yline(7.3523, lpattern(dash)) connect(direct) recast(line) ytitle("LR Statistics") xtitle("First Threshold") name(LR21)
_matplot e(LR22), columns(1 2) yline(7.3523, lpattern(dash)) connect(direct) recast(line) ytitle("LR Statistics") xtitle("2nd Threshold Parameter") name(LR22)
graph combine LR21 LR22, cols(1)